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Anisotropic emission of gravitational-waves (GWs) from inspiralling compact binaries leads to 
the loss of linear momentum and hence gravitational recoil of the system. The loss rate of linear 
momentum in the far-zone of the source (a nonspinning binary system of black holes in quasi- 
circular orbit) is investigated at the 2.5 post-Newtonian (PN) order and used to provide an analytical 
expression in harmonic coordinates for the 2.5PN accurate recoil velocity of the binary accumulated 
in the inspiral phase. The maximum recoil velocity of the binary system at the end of its inspiral 
phase (i.e at the innermost stable circular orbit (ISCO)) estimated by the 2.5PN formula is of the 
order of 4kms _1 which is smaller than the 2PN estimate of 22kms _1 and indicates the importance 
of higher order post-Newtonian (PN) corrections. Going beyond inspiral, we also provide an estimate 
of the more important contribution to the recoil velocity from the plunge phase. The maximum 
recoil velocity at the end of the plunge, involving contributions both from inspiral and plunge phase, 
for a binary with symmetric mass ratio v = 0.2 is of the order of 182 km s -1 . 

PACS numbers: 04.25. Nx, 04.30.Db, 97.60.Jd, 97.60.Lf 



I. INTRODUCTION 

A coalescing black hole (BH) binary which is anisotropically emitting gravitational-waves (GWs) will experience a 
recoil as a consequence of the loss of linear momentum from the binary through outgoing GWs. This phenomenon 
of gravitational- wave recoil has substantial importance in astrophysics especially if one wants to study models which 
suggest the formation and growth of super massive black holes (SMBHs) at the centers of galaxies through successive 
mergers from smaller BHs (stellar or intermediate mass BHs) PJ . If the kick velocity of the product BH is more than 
its escape velocity from the host galaxy, the formation of SMBHs will not be favored, as would be the case with dwarf 
galaxies and globular clusters (see Q for observational evidence for ejection of the SMBHs). Even if the recoil velocity 
of the product BH is not sufficient to eject it from the host (which may be the case with giant elliptical galaxies), the 
product BH would be displaced from the center and eventually would fall back. Such a process may have important 
dynamical changes at the galactic core. For a more detailed overview of astrophysical possibilities, see Ref.M. One 
important claim of [l| is that models that grow SMBHs from the mergers of smaller BHs will not be favored for the 
galaxies at redshifts z > 10 due to the difficulty in retaining the 'kicked' black holes. However, observations of the 
local universe suggest that most of the galaxies (more than 50% of them) have SMBHs at their centers [||. This shows 
that there must be something which prevents the ejection of the central black hole from many of these galaxies and 
Ref. Q investigates such questions. In the light of the above arguments it becomes important to have an accurate 
estimate of the recoil velocity of the coalescing binary black holes (BBHs). 

One of the first proposals investigating the phenomenon of gravitational- wave recoil is due to Peres @. Fol- 
lowing the analogy from the classical electrodynamics it was suggested that the lowest order secular effects 
related to gravitational-wave recoil arise due to the interaction of mass quadrupole moment with mass octupolc 
moment or current quadrupole moment. His work provided the first formal theory for gravitational- wave recoil 
of a general material system in linearized gravity and is valid for any kind of motion (rotational, vibrational 
or any other kind) given the source is localized within a finite volume. In another early work, Bonnor and 
Rotenberg [|| studied the emission of gravitational-waves from a pair of oscillating particles and suggested the 
possibility of GW recoil. Papapetrou [7j derived the leading order formula which involved interaction of mass 
quadrupole moment with mass octupole moment and current quadrupole moment. Later, Thorne Q generalized 
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the idea by providing a general multipole expansion for the linear momentum loss as seen in the far zone of the source. 

Within the post-Newtonian (PN) scheme, the leading order contribution to the linear momentum flux from an 
inspiralling binary system of two point masses in Keplarian orbit was computed by Fitchett @ and binary motion in 
circular orbit was discussed as a limiting case of the main results of the work. The first PN correction was added to it 
by Wiseman and the circular orbit case was discussed as a special case. In a work by Blanchet, Qusailah and Will 
(hereafter BQW) (TT| the PN expression for linear momentum flux, from a nonspinning compact binary moving in 
quasi-circular orbit, was extended to the 2PN order by ad ding hereditary contribution that occurs at 1.5PN order and 
the instantaneous one occurring at the 2PN order. Kidder [12| computed the leading order (spin-orbit) contribution to 
linear momentum flux for generic orbits and discussed circular orbit effects as a limiting case of his findings. Recently, 
Racine et al. [l3| extended Kidder's work by adding higher order spin corrections (spin-orbit terms at 1.5PN order, 
spin-orbit tail and spin-spin term at 2PN order). They provided 2PN accurate expression for the linear momentum 
flux from spinning BBHs in generic orbits. They also specialize to the binary motion in circular orbit and provide 
estimates for the recoil velocity, accumulated during the inspiral phase, for equal mass binaries with spins equal in 
magnitude but opposite in direction. 

Using the black hole perturbation (BHP) theory, Favata et al. [l4[ estimated the recoil velocity of the binary, 
treating it as a test particle inspiralling into a black hole (spinning or nonspinning), up to the inner- most stable 
circular orbit (ISCO) (accounting for the recoil velocity accumulated during the inspiral phase) very accurately. 
Though their calculations were valid only in extreme mass ratio limit (q m = m\/m2 << 1) they extrapolated their 
results to g m ~ 0.4 (with modest accuracy) using some scaling results from quadrupole approximations. A crude 
estimate of the contributions due to the plunge was also given. Within the validity of the approach, their estimates 
suggested the typical recoil velocity can be of the order of 10-100 kms -1 but for some configurations it may reach 
rou ghly up to 500 kms -1 . Another computation by Damour and Gopakumar (l5j used effective one-body approach 
0,113 to compute the total recoil velocity of the final black hole taking into account the contributions from all the 
three phases (inspiral, plunge and ringdown). Depending upon the method they used to compute linear momentum 
flux their estimates for maximum recoil velocity lie in the range 49-172 kms -1 . Ref. [l8| presents estimates of the 
recoil velocities for binaries in orbits with small eccentricities using an approximation technique that is valid only for 
late stages of the plunge. They also combine their results with the PN estimates of recoil velocity at ISCO of BQW in 
order to give estimates for recoil velocity for binaries in quasi-circular orbits and find that for a binary with symmetric 
mass ratio (ratio of reduced mass of the binary to the total mass) v — 0.2 the recoil velocity estimates should lie in the 
range of (79-216) kms -1 . In a recent work [l9[ the recoil of the final BH was investigated combining the results of [ll[ 
with the calculation of contribution from the ringdown phase performed using close-limit approximation. They found 
that the radiation emitted in the ringdown phase produces a significant anti-kick and thus brings down the estimates 
of recoil velocity based on only inspiral and plunge phase e.g. after including the contributions from ringdown phase 
the maximum recoil velocity of the final black-hole is of the order of 180 kms -1 as compared to BQW estimate of 
243 kms -1 which does not include the contribution from the ringdown phase (also see Fig.l of [l^| for a comparison 
of this result with various numerical and analytical estimates) . In another recent work [20j the phenomenon of recoil 
of a spinning BBH (extreme mass ratio) due to the inspiral, merger and ringdown phase of its evolution has been 
investigated. The issue of antikick has been examined very carefully and they found that for orbits aligned with the 
BH spin, the antikick grows with the spin. Also, a prograde coalescence of a smaller BH into the rapidly rotating 
bigger BH results in the smallest kick whereas the retrograde coalescence insures the maximum recoil. 

In addition to the analytical or semi-analytical estimates of the recoil there have been many investigations using 
numerical techniques. Recent numerical simulations for nonspinning (2l|-|24| BBHs in quasi-circular orbit have shown 
that the recoil velocity can reach up to a few hundred kms -1 while for spinning case (25l - [28j the recoil velocity 
estimates are much higher and can be of the order of few thousand kms -1 . Although numerical simulations can put 
better constraints on these estimates, such simulations (especially those which includes BH spins) are computationally 
very expensive. Moreover a very detailed multipolar study of numerical results for BBH recoil [29| shows the need of 
analytical and semi-analytical schemes in order to gain a deeper understanding of the problem at hand and also as a 
check to numerical results. 

In the present work we extend the 2PN calculation of fill ] for linear momentum loss from a nonspinning BBH in 
quasi-circular orbit by adding terms (both instantaneous and hereditary) which contribute at 2.5PN order and thus 
give an analytical expression for linear momentum flux which is now 2.5PN accurate. Naturally in the 2PN limit our 
expression for linear momentum flux given by Eq. (|3.23|) reduces to Eq. (20) of [111. The 2.5PN accurate expression for 
the recoil velocity of the binary is given by Eq. (|4.6[) which reduces to Eq. (23) of [Tl| in the 2PN limit. For computing 
the contribution to the recoil velocity due to the plunge phase we simply adopt the discussion given in Sec. (4.1) 
of [ll] ] and perform the computation using our 2.5PN accurate formulae. We find that the maximum recoil velocity 
experienced by the binary at the end of inspiral (at ISCO) and end of the plunge (which includes the contributions 
from both inspiral and plunge phase) is of the order of 4 kms -1 and 182 kms -1 respectively and corresponds to the 
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binary with symmetric mass ratio of v = 0.2. In contrast, the maximum recoil velocity at the end of the inspiral 
and the plunge using the 2PN formulae (see Fig. 1 of |Tl|) is of the order of 22kms _1 and 243 km s -1 corresponding 
to the same v = 0.2. We see here that inclusion of terms at 2.5PN order brings down the estimates for the recoil 
velocity significantly, exhibiting in this problem the feature arising from the asymptotic nature of the PN expansion 
and the need to exp licitly investigate the next PN order implications of a calculation. This also reminds us of a 
similar result of [H| where the inclusion of 1PN contribution brought down the Newtonian estimates since the 1PN 
term contributed negatively to the recoil velocity. Something similar happens here and the large negative coefficients 
at 2.5PN order (see Eq. (|4.6|0 brings down the 2PN estimates significantly. 

The paper is organized in the following way. In Sec. [TT| we hrst provide the PN structure of the linear momentum 
flux in terms of the radiative multipole moments and then we give explicit expressions for the instantaneous and 
hereditary contribution separately in terms of the source multipole moments. Section IlIII starts with the formulae 
for source multipole moments with desired PN accuracy and next shows the computation of both instantaneous and 
hereditary contributions to the linear momentum in the far-zone of the binary. Finally, we give the 2.5PN accurate 
expression for the linear momentum flux by adding instantaneous and hereditary contributions. In Sec.|lV]we discuss 
the computation of the recoil velocity of the binary and also give the 2.5PN accurate analytical expression for the 
same. Section [V] explores the method for estimating the recoil velocity accumulated during the plunge phase. In 
Sec. IVII we present our numerical estimates of total recoil velocity and its dependence on the composition of the 
binary as well as final discussions. 



II. THE POST-NEWTONIAN STRUCTURE FOR LINEAR MOMENTUM FLUX 

The general formula for linear momentum flux in the far-zone of the source in terms of symmetric trace-free (STF) 
radiative multipole moments is given in Q and at relative 2.5PN order it takes the following form (see e.g. Eq. (4.20') 
of Ref.@].) 
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In the above expression Ul and Vl (where L = ■ ■ ■ ii represents a multi-index composed of I spatial indices) are 



the mass-type and current-type radiative multipole moments respectively and Ut^ and V} p ' denote their p th time 
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derivatives. The Levi-Civita tensor is denoted by e^fc, such that 



-123 — +1 and 0(l/c 6 ) indicates that we ignore 
contributions of the order 3PN and higher. The moments appearing in the formula are functions of retarded time 
U = T — (R/c) in radiative coordinates. Here R and T denote the distance of the source from the observer and the 
time of observation in radiative coordinates, respectively. 

It should be evident from Eq. (|2.1j) that the computation of 2.5PN accurate linear momentum flux requires the 
knowledge of Uij,Vij and Ujjk at 2.5PN order, Ujki and Vijk at 1.5PN order, and Vijkim and Vijki at Newtonian 
order. In a recent work [301 ] . Ul and Vl have been computed with accuracies sufficient for the present purpose 
using multipolar post-Minkowskian (MPM) approximation approach ITf 36 . In the MPM formalism Ul and Vl arc 
related to canonical moments Ml and Sl (Eqs. (5.4)-(5.8) of [3(3]) which in turn are related to source moments 
{II-, Jl, Xl,Wl,Yl, Zl} (Eqs. (5.9)-(5.11) of [30]). Rewriting the expressions for the radiative moments in terms of 
the source moments, the linear momentum flux can be decomposed as the sum of two distinct parts: the instantaneous 
terms and the hereditary terms. By instantaneous we refer to contributions in the linear momentum flux which depend 
on the dynamics of the system only at the retarded instant U = T — (R/c). Hereditary contributions to the flux, on 
the other hand, are terms nonlocal in time depending on the dynamics of the system in its entire past 37J. The linear 
momentum flux thus is conveniently decomposed into 
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where the instantaneous part is given by 
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and the hereditary contribution can be written as 
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Here M denotes the ADM mass of the system. To appearing in above hereditary integrals is an arbitrary constant 
and and is related to an arbitrary length scale, ro, by the relation To=ro/c. It enters the relation connecting retarded 
time U = T — (R/c) in radiative coordinates to retarded time u = th — r\ x jc in harmonic coordinates (where rh is the 
distance of the source in harmonic coordinates). The relation between retarded time in radiative coordinates, and the 
one in harmonic coordinates reads as 
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III. THE 2.5PN LINEAR MOMENTUM FLUX: APPLICATION TO INSPIRALLING COMPACT 

BINARIES IN CIRCULAR ORBITS 



Equations ()2.2p - (|2.5j) collectively give the far-zone linear momentum flux from generic PN sources in terms of 
the source multipole moments {II, Jl, Wl, Yl}. In this section we specialize to the case of nonspinning inspiralling 
compact binaries, in which two compact objects (neutron stars and/or black holes) are moving around each other in 
quasi-circular orbits. All the source multipole moments in case of nonspinning inspiralling compact binaries moving 
in quasi-circular orbits are now known with the accuracies sufficient for the present purpose and have been computed 
and listed in [3(j (see Eqs. (5.12)-(5.25) there). Here we just quote those results with the accuracies that is required 
for the present purpose. For mass-type moments we have 
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and, for the current-type moments 1 we have 
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Computation of linear momentum flux at 2.5PN order also requires, 1PN accurate expression for mass monopolc 
moment, which can be identified with ADM mass (M) of the source, and Newtonian accurate expression for the 
current dipole moment Jj. We have 
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In addition to mass-type and current-type moments we also need some of the gauge moments which only need to be 
Newtonian accurate and are given as 
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In the above, m = mi +TO2 is the total mass of the binary with mi and 777,2 as binary's component masses and v is the 
symmetric mass ratio and is defined by the combination (mi m2/m 2 ). x 1 and v l denote binary's relative separation 
and relative velocity of the two objects constituting the binary respectively and can be defined as x % = y\ — y\ and 
v % = Ax 1 1 At = v\ — v\ (where (2/1,2/2) an( i ( v \i v 2) are center-of-mass positions and velocities of components of the 
binary). Angular brackets, (), surrounding indices denote symmetric trace- free projections. 7 is a PN parameter and 
is defined by the quantity (Gm/c 2 r). 



A. Instantaneous Terms 



Equation (|2.3[) is the general formula for the instantaneous part of the linear momentum flux from generic PN 
sources in terms of the source multipole moments {II, Jl, Wl, Yl}- Computation of linear momentum flux involves 
computing time derivatives of the source multipole moments which in turn requires the knowledge of equations of 
motion with appropriate PN accuracy. Linear momentum flux computation at 2.5PN order will thus require 2.5PN 
accurate equations of motion (30l[38|. 

Let x-y plane be the orbital plane of the binary. 2 <fi being the orbital phase of the binary giving the direction of 
the unit vector, n = x/r, along the binary's relative separation by 

n = cos <j> e x + sin cj> e y . (3-5) 



1 The coefficient "-484/105" appearing at the 2.5PN order in the expression for Jij in Eq. 5.15b of l3dH is incorrect and should be replaced 
by "-148/35" (see JOijJ). 

2 Since we are considering only nonspinning binary systems in quasi-circular orbits, the motion will be in a fixed plane. 
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The binary's relative separation, velocity and acceleration are given by 

x = r n, (3.6a) 
v = fn + rwA, (3.6b) 
a = (r — r a; 2 ) n + (rii) + 2fw) A. (3.6c) 

where an over dot denotes a time derivative and r = |x| is the distance between the two objects in the binary. The 
orbital frequency lo is given by uj — <f>. The motion of the binary can be described by the rotating orthonormal triad 
(n, A, e z ) with A = e z x n. 

Up to 2PN order one can model the binary's orbit as exact circular orbit with r = 6j = 0, but at 2.5PN order 
orbit of the binary decays due to radiation reaction forces and one must include the inspiral effects. The leading 
order effect is computed using energy balance equation assuming that system is losing its orbital energy only through 
gravitational radiation. At the 2.5PN order, for f and <j> we have 
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By substituting the expressions for r and Cj in Eq. (j3.6bp ~ p.6cp one can write for the relative inspiral velocity and 
relative acceleration as 
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Finally we give the PN expression for orbital frequency as a function of the binary's separation r which is now known 
with 3PN accuracy (30j but in the present work we just need the 2PN accurate expression. In harmonic coordinates 
it is given as 
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It is often convenient to use a parameter x which is directly connected to the orbital frequency rather instead of using 
the PN parameter 7. 3 Our new parameter x is related to orbital frequency (Eq. 6.5 in [30| ) as 

n \ 2 / 3 



(3.10) 

A relation between 7 and a; can be obtained by using Eq. (|3.10p in Eq. (|3.9p and inverting for 7 in terms of x. At the 
2PN order, the PN parameter 7 is related to the parameter x as 

7 = .{l + x(l-9+, 2 (l-^)+o( T i)}. (3.11) 

Now we have all the inputs to compute the derivatives of source multipole moments with accuracies sufficient for the 
computation of 2.5PN accurate expression for linear momentum flux. Once we have the desired time derivatives of 
various source multipole moments we can insert them in Eq. (|2.3p to get the 2.5PN accurate instantaneous part of 
the linear momentum flux. After a tedious but straightforward computation we get for the 2.5PN accurate expression 
for linear momentum flux in terms of the parameter 7 
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The use of a; as a PN parameter is useful since it remains invariant for a large class of coordinate transformations including the harmonic 
and Arnowitt, Deser and Misner (ADM) coordinate systems. 
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Alternatively, we can rewrite the instantaneous part of linear momentum flux given by Eq. (|3.12l) in terms of the 
parameter x by using Eq. (|3.11[) in the above equation. At 2.5PN order the linear momentum flux in terms of the 
parameter x reads as 
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B. Hereditary Terms 

In this subsection, we shall compute the hereditary contribution to linear momentum flux from a nonspinning 
inspiralling compact binary in quasi-circular orbits, which in terms of the source multipole moments is given by 
Eq. (|2.5[) . The leading order hereditary contribution (1.5PN term) for the nonspinning compact binaries in circular 
orbit have been computed in [ll[ and later confirmed by Racine et al. [l3| . In this section we extend the computation 
of the hereditary contributions by adding terms contributing at the 2.5PN order. 

If x-y plane is the binary's orbital plane and the orbital phase at a given retarded time U be (f>(U) then unit vectors 
n and A can be written as 



n(U) = cos (j)(U) e x + sin <j)(U) e y , 
X{U) = —sin (f>(U) e x + cos 4>{U) e y 



(3.14a) 
(3.14b) 



It is evident from Eq. (|2.5[) that to compute the hereditary contribution one must know the relevant multipole moments 
and their derivatives both at any retarded time U as well as at some other time U' = U — r < U. Since multipole 
moments at retarded time U' shall involve n and A at U', it would be useful to express n(t/') and X(U') in terms of 
n(/7) and X(U), which are independent of the integration variable, r, and thus one can pull out the vector quantities 
out side the hereditary integral. Following [l3|, one possible way is to express h(U') and X(U') as a linear combination 
of n(U) and X(U) as 



h(U') = cos [<j)(U) 
X(U') = sin [4>(U) - 



■4>(U')]n-sm [</>(U)-<f>(U')]\, 
<f>(U')] n + cos [<t>{U) - <j)(U')] A . 



(3.15a) 
(3.15b) 



It should be evident from Eq. (|2.5|) that hereditary contributions at the 2.5PN order require 1PN accuracy for the 
quantities appearing in first four terms while in remaining six terms they need only be Newtonian accurate. It should 
be clear that while computing the time derivatives of the source multipole moments for hereditary contributions, the 
equations of motion need be only 1PN accurate at most. To start, let us consider the combinations of derivatives of 
source multipole moments appearing in the first term of Eq. (|2.5[) . i.e l\fl{U) 1^' \U — t). We can use Eq. (|3.15l) to 

express the quantity at hand in terms of n(U) and X(U) and then perform the contraction of indices. After some 
straight forward algebra we have 
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where we have defined 5(f) = (f>(U) — <f>(U — r). 

Similarly we can write for combinations of source multipole moments in various terms of Eq. (|2.5[) as 
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(3.17a) 
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(3.17h) 
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It is evident from the above that the dependence of the relevant quantities on the integration variable r is only 
through 50 which under the assumption of adiabatic inspiral takes the form 



0(£7) - 0((7 - r) 
0(10 



T = U 



where second and higher derivatives of have been neglected. 

Finally, one just needs the following standard integral to compute the hereditary terms in p. 5 



log 
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(3.18) 



(3.19) 



Eqs. (|3.16|) - (|3.19[) provide all the necessary inputs that are needed for computing the hereditary terms. For the 
sake of compactness of the paper we wish to skip some of the intermediate outcomes of the calculation and directly 
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quote the final expression for the 2.5PN accurate hereditary contribution which in terms of the parameter x reads as 
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(3.20) 



where ujq appearing in the above provides a scale to the logarithms and is given as 



1 /5921 48, n 405, „ „ 
— exp 1 log 2 log 3 — C 
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(3.21) 



where C is Euler's constant. One can verify that, terms involving the logarithms of frequency log (^j^j appearing in 

Eq. (|3.20p can be reabsorbed into a new definition of phase variable and thus will disappear from the final expression 
for linear momentum flux. This possibility of introducing a new phase variable containing all the logarithms of 
frequency has been noticed and used in earlier works [TTI . |39i . liTj ] . We define the new phase variable tp as 



, 2GMuj , / uj 

lp = 4> 3 log 

C 6 \LOqi 

where M is the ADM mass of the source and is given by Eq. (|3.3ap . 



(3.22) 



C. Total LMF 



The final expression for the LMF in terms of the parameter x can be obtained by simply adding Eq. (|3.13[) and 
Eq. (|3.20p and using ip, given by Eq. (|3.22p . as the phase variable. In the final form the 2.5PN expression for LMF 
reads as 
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(3.23) 



V 12615 841 ° 841 
It should be clear that now n and A are in the direction of new phase angle tp and tp + n/2 respectively and are given 



n = cos ip e x + sin ip e y , 
A = —sin ip e x + cos ip e y . 



(3.24a) 
(3.24b) 



where ip is given by Eq. (l3.22l) . 



IV. RECOIL VELOCITY 



Given the 2.5PN far-zone linear momentum flux due to a nonspinning inspiralling compact binary in quasi-circular 
orbits (Eq. (|3.23p ) one can have 2.5PN accurate formula for the loss rate of linear momentum by the source using the 
linear momentum balance equation, which is 
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The net loss of linear momentum can be obtained by integrating the balance equation, i.e. 

t-i 



AP l 



dt' r P 



Since we have for the nonspinning compact objects moving in quasi-circular orbit 

dh l u 

-dt = ~" n > 



where ui is the orbital frequency of the inspiral. We can have 

h l dt -- 
Xdt -- 
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— dn=— n 

UJ UJ 



(4.2) 

(4.3a) 
(4.3b) 

(4.4a) 
(4.4b) 



It should now be clear that in order to compute the net loss of linear momentum at any instant (before the orbit 
evolves to the last stable orbit (LSO) around which generally the assumption of adiabatic approximation should break 
down) one just needs to replace n l by —\/uj and X 1 by h l /uj. Once we have the net change in the momentum during 
the orbital evolution of the binary we can obtain the recoil velocity of the source by simply dividing it by the mass of 
the system i.e. 



AV l = APVm , 

and we find in terms of our parameter x, the 2.5PN accurate expression for the recoil velocity as 



(4.5) 
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(4.6) 



NUMERICAL ESTIMATES OF RECOIL VELOCITY 



With 2.5PN accurate formulae for the linear momentum flux (Eq. Q3.23p ) and the recoil velocity (Eq. (|4.6|) ) we would 
now wish to compute the recoil velocity accumulated during the plunge phase. Generally, the PN approximation is 
considered to be less reliable for the orbits within the ISCO, by this we mean that PN corrections, when compared 
to the leading order contribution, become comparable. If these corrections are small even beyond the ISCO then one 
can use Eq. (|3.23j) to estimate the velocity accumulated during the plunge phase. Generally it is expected that these 
corrections would become comparable to the leading order contribution near the common event horizon and thus one 
can only provide a crude estimate of the recoil velocity accumulated during the plunge. For this purpose we simply 
adopt the methodology used in BQW [ll[ (see Sec. 4.1 there 4 ). We shall first compute the recoil velocity at the ISCO 
using Eq. (|4.6[) , where ISCO is taken to be that of a point particle moving around a Schwarzschild black hole with 
the mass equal to the total mass of the binary i.e. m = m\ + m?,. For the kick velocity at the ISCO we can write 



Vr 



464 



ISCO 



105 
147101 
68904 



c V 1 — Av x 



ISCO 



5/2 
C ISCO 



1 + XiscO 



2663 
116 



2185 
~87~ 



452 
'~87 



■ 7T V 



1139 
522 



309 3/2 
58 7rXlsco 



^ISCO 



71345 36761 



22968 2088 



5/2 
C ISCO 



106187 
50460 



32835, „ 77625, n 
■ log 2 + log 3 



841 



3364 



4 Though all necessary details have been given in [lll | we provide here some of the basic formulae for completeness of the paper as well 
as for the convenience of the reader. 
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where xisco 5 "isco ano - ^isco d en °te values of x,h l and A 2 at the ISCO respectively. For a point mass moving 
around a Schwarzschild black- hole of mass m in circular orbits, xisco=l/6. We choose for simplicity the phase at 
the ISCO to be ip = and thus n\ sco = {1, 0, 0} and Aisco = {0, 1, 0}. With this we can compute the recoil velocity 
due to the inspiral phase. 

Following BQW, we adopt the effective one-body approach [T^, [T3] to compute the plunge contribution to the recoil 
velocity. We assume that a point particle of mass fx is moving in the gravitational field of a Schwarzschild black hole 
of mass m where n is the reduced mass of the system. In addition to this we also ignore the effect of radiated energy 
and angular momentum on the plunge orbit. We have for the geodesic equations in the Schwarzschild geometry as 

(5.2a) 
(5.2b) 
(5.2c) 



dt 
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dr s \ ,=,9,9 9 / „ 2Gm\ I L 



dr ) \ c 2 r s 

Here r is the proper time along the geodesic and E and L are energy and orbital angular momentum per unit mass 
and can be defined in terms of dimensionless variables E and L as 

E = c 2 E, (5.3a) 
L=—L. (5.3b) 



c 

Using Eq. (|5.2bp - (|5.2c[) one can obtain the phase of the orbit as 



I" l T— —-7TT 2 \ dy, (5.4) 



!I0 



[£ 2 -(l-2 2,)(l + LV)] 



with y = (Gm/r s c 2 ) and for the phase at the beginning of the plunge we have chosen, ip = (at y = yo), to match 
the phase of the orbit at the ISCO. 

Now the kick velocity accumulated during the plunge phase can be given by the following formula 

I /"^Horizon Jpi 

where to and tHorizon are the times at the beginning of the plunge and when particle approaches the horizon, respec- 
tively. 

Now it is clear from Eq. (|5.2a|) that the time coordinate t is singular in nature at the horizon i.e. at r s = 2Gm/c 2 
and thus we must have a variable which is nonsingular in nature at the horizon to compute the above integral. We 
can write Eq. (|5.5|) as 

1 r ( du) \ dP 

'plunge - — J [d^Jdi J ~dt ' 

where the quantity, Q, is the proper angular frequency, defined asw = dip/dr, and is given by Eq. (|5.2b[) . After some 
straightforward algebra wc obtain 

_ Gl_ rv**~n f l dP l \ ( dy \ 

^W" C 3 J m [ x S/2 dt ) { [ E2_ {1 _ 2y){1 + L 2y2 ) f 2 ) ' ^ 

where dP/dt in terms of our parameter x is given by Eq. (|4.1|) in combination with Eq. ()3.23j) . Now x is related to 
the variable y by the following equation and can be obtained after few steps of algebra 



AK\ ungc = - — — — (5.6) 



f "1 2/3 

| y 2 (l-2y) 



(5.8) 
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Using the above relation and the definition of phase given by Eq. (|5.4|) , the quantity inside the integral of Eq. (|5.7|) 
becomes a function of just the integration variable, y. With known values of E and L for the plunge one can 
numerically compute the integral of Eq. ()5.7j) . 

Our next task is to choose appropriate values for E and L which are also consistent with their values at the ISCO. 
In fact there may be several ways to match a circular orbit at the ISCO to a suitable plunge orbit; we would use the 
two methods which have been used in [ill ]. In the first method, the particle is given an energy E = c 2 E such that, 
at the ISCO, and for an ISCO angular momentum Zigco = %/l2 (Gm/c), its radial velocity is given by the standard 
quadrupole energy-loss formula for a circular orbit, which is given as 



dr h _ 64 /Gm x 
dt 5 V c 2 rh 



vc, (5.9) 



where i\ is the binary's radial separation in harmonic coordinates. For a test particle, at the ISCO, = 5 (Gm/c 2 ), 
so we have (dr^/ 'dt)isco = — (8/25) 2 v c. Since radial and time coordinates in Schwarzschild and harmonic coordinate 
systems are related as 

r s = r h + G m/c 2 , t s = f h = t . (5.10) 
we have for the radial velocity of the particle in Schwarzschild coordinates as 

G?r s dry. 

— - = — - 5.11) 

dt dt K ' 

= - (8/25) 2 vc. 

It is easy to show using Eq. (|5.2aj) and (|5.2cj) that the required energy for such an orbit will be given by 
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where dr s /dt is given by Eq. (|5.11[) . 



Now with E given by above equation and the choice L = vl2 we can compute the desired integral numerically. 5 
As last input, for the limiting values of the integration variable, y we choose, yo = 1/6 and ?/Horizon = (2(1 + v)) . 
Note that the choice yo=l/6 and the one that have been made for E and L above will not be consistent with Eq. (|5.8[) 
and thus when computing the recoil velocity at the ISCO the value of the parameters x at ISCO must be consistent 
with the choice for yo, E and L made above. 

In the second method, one matches the circular orbit at the ISCO and the one associated with the plunge by 
evolving it across the ISCO and it can be performed using the energy and angular momentum balance equations for 
circular orbits in the adiabatic limit at the ISCO. For this we shall have 

dE 32 c-v 5 

■x ISCO , (5.13aJ 



dt 5 G m 

dL ( Gm wisccA 1 dE 32 c 3 v 7/2 
~dt = { c 3 ) ~dt=~T~G^ Xlsc °- 



(5.13b) 



: m 

Following [llj , we can write for the quantity on the left side of Eq. (|5.13[) 

dE/dt = (E — E ISCO )/(a P), (5.14a) 

dl/dt = (L- L lsco )/(a P). (5.14b) 

Here a denotes a fraction of the orbital period P of the circular motion at the ISCO. Now using, wisco = 
(c 3 /Gm) aifgQo we haw for the plunge orbit, 



E = Eisco - -^-a^^isco > ( 5 - 15 ) 



5 For this purpose we shall use the NIntegrate option inbuilt in Mathematica. 
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FIG. 1: Recoil velocity as a function of the mass parameter v (symmetric mass ratio) has been shown. Plot shows a comparison 
between the results obtained using the two methods (we call them Ml and M2), discussed in Sec|V] that have been used to 
match the circular orbit at the ISCO to a suitable plunge orbit. It is evident from the figure that the two methods are consistent 
with each other for smaller values of the parameter a. 
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Finally, in the second model, in order to integrate the integral in the problem we need to specify the limiting values 
for the variable y. For the initial value of the parameter (at ISCO), y — yo, one can solve the following equations 
which is obtained using Eq. 



a^isco = 6~ 



2/o(l-2y ) 



2/3 



(5.17) 



For the value of the parameter, y, at the horizon again we take, ^/Horizon = (2(1 + v)) . For the fraction a of the 
period, we choose values between 1 and 0.01, and check the dependence of the result on this choice (see below). 



VI. RESULTS AND DISCUSSIONS 



Equation (|3.23[) gives the 2.5PN formula for computing the loss rate of linear momentum in the far-zone of a 
nonspinning inspiralling compact binary in a quasi-circular orbit. In Scc|V]we show how one can numerically estimate 
the recoil velocity accumulated during the plunge phase, after making some simplifying assumptions. The recoil 
velocity at the end of the inspiral phase, i.e. at the ISCO, is given by Eq. (|5.1[) whereas the recoil velocity accumulated 
during the plunge phase is given by Eq. (|5.7p . The integral of Eq. (|5.7|) needs to be evaluated numerically keeping 
in mind that appropriate choices for energy and angular momentum at the onset of the plunge phase has been made 
in order to match the inspiral and plunge orbits at the ISCO. Figure Q] shows our numerical estimates for the recoil 
velocity, based on the two methods (we call them Ml and M2) for matching the circular orbit at the ISCO to a 
suitable plunge orbit, discussed in the previous section. Figure. [T] also shows a comparison between the recoil velocity 
estimates using the two methods, Ml and M2. It is evident from the figure that results from both the methods are 
consistent with each other for smaller value of the parameter a, defined above. In case of M2, we have shown curves 
corresponding to a = {0.01,0.05,0.1,0.5,1.0} and one can see that the curves with a = 0.01and0.05 are very close 
to the curve corresponding to the Ml. We also observe that the maximum recoil velocity, as shown in Fig. [T] is 
~182kms~ 1 (for Ml and M2 with a = 0.01) and corresponds to the value of the parameter v=Q.2. This is lower 
than the 2PN accurate BQW estimate of about 243 km s -1 for the binary with the same mass ratio (i.e. v—0.2). This 
behavior is due to the presence of large negative coefficients at the 2.5PN order (see Eq. (|4.6[) ) which bring down the 



15 



estimates significantly. Such a behavior is not new to PN calculations, e.g. a similar behavior was observed in [l(| 
at 1PN order (see Fig.l of [H|), where the use of 1PN accurate results give a lower estimate for the recoil velocity as 
compared to the one obtained using the Newtonian formulae since the 1PN term again contributes negatively to the 
recoil velocity. 

Note that the estimates presented in Fig. [1] use only the leading order radiation reaction effects for setting initial 
energy (in Ml) and energy and angular momentum (in M2). We repeat the exercise using 2.5PN expressions for 
relevant quantities beyond the leading order effect and find that changes in estimates are negligible (relative % 
changes are less than 0.5%). 

As discussed earlier, normally the PN approximations are expected to become less and less reliable beyond the 
ISCO. This leads to a crude estimate of the accumulated recoil velocity during the plunge phase. Hence, it becomes 
important to compare our results to some other numerical/analytical estimates, in order to be sure that these estimates 
are indeed reliable. In the case of the present work the closest comparison for the recoil velocity estimates can be 
made by comparing our results with those of BQW [ll|. For a binary with v = 0.22 and v = 0.23, BQW suggest 
that the recoil velocity should lie in a range, (171-251) kms -1 and (146-220) kms -1 , respectively. The uncertainty 
in their results has been estimated by flexing the 2PN expressions by addition of 2.5PN, 3PN and 3.5PN terms and 
then computing the maximum variation in their results (see llj for details). Our estimates of the recoil velocity for 
a binary with v = 0.22 and v = 0.23 are 171 kms -1 and 158 kms -1 , respectively, and thus our estimates lie in the 
window for the recoil velocity provided by BQW. However, we should note here that our estimates can also change 
if we add contributions coming from the 3PN and the 3.5PN terms (although changes may be relatively smaller). 
Currently, such an extension is not possible as we do not have sufficiently accurate inputs in order to perform such 
computations and thus it will be the subject matter of a work in the future. Our estimates are also consistent with 
an earlier numerical work [2l[ which suggests a range of values for recoil velocity between (100 — 380) km s -1 and 
(90 — 290) kms -1 for v = 0.22 and v = 0.23, respectively. As discussed in Sec. HI Ref. [IH suggests that maximum 
recoil velocity estimate for a binary with v = 0.2 in quasi-circular orbit lie in a range between (79-216) kms -1 . As 
mentioned above, our estimate for such system is 182 kms -1 and thus is consistent with their estimates. 

We witnessed above that inclusion of 2.5PN contributions significantly changed earlier PN estimates for the recoil 
velocity indicating that contributions at higher orders need to be explicitly assessed due to the asymptotic nature 
of the PN expansion. As mentioned above, contributions at other high PN orders such as at 3PN and 3.5PN order 
should be included in some future work in order to have better estimates for the recoil velocity, although the changes 
may be relatively smaller as compared to those brought in by 2.5PN contributions. A numerical study [22[ suggests 
that the recoil velocity estimates at the ISCO should be of the order of ~ 14 kms -1 for a binary with v = 0.24 and this 
estimate matches well with BQW estimates for the same system. This is a relatively higher estimate as compared to 
our estimate of 2.2 kms -1 at the ISCO for a system with the same mass ratio. In such a case, we should expect that 
inclusion of higher order contributions at the 3PN order will contribute to the recoil velocity positively (in contrast 
to the negative contributions form 2.5PN terms) and thus could bring up the estimates to match with estimates of 
[22l | and BQW. In addition to this, as a follow up of this work, one can try to include contributions due to the final 
ringdown phase using 2.5PN accurate initial conditions 6 and then combine this with the recoil velocity estimates for 
the inspiral and plunge phase presented here. This will allow one to make more direct comparisons with the results 
obtained using numerical relativity and effective one-body approach which include contributions from all three phases 
of the binary evolution. 
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